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Abstract 



An ab initio Wannier-function-based approach to electronic ground- 
state calculations for crystalline solids is outlined. In the framework of 
the linear combination of atomic orbitals method the infinite character 
of the solid is rigorously taken into account. The Hartree-Fock ground- 
state energy, cohesive energy, lattice constant and bulk modulus are 
calculated in a fully ab initio manner as it is demonstrated for sodium 
chloride, NaCl, using basis sets close to the Hartree-Fock limit. It 
is demonstrated that the Hartree-Fock band-structure can easily be 
recovered with the current approach and agrees with the one obtained 
from a more conventional Bloch-orbital-based calculation. It is argued 
that the advantage of the present approach lies in its capability to 
include electron correlation effects for crystalline insulators by means 
of well-established quantum chemical procedures. 

1 Introduction 

Calculations of ground-state properties as well as band-structure calcula- 
tions are usually done in the so-called Bloch-orbital approach that employs 
itinerant orbitals, which are deiocalized all over the infinite crystal]]]]. The 
wavefunctions obtained, in a single-determinantal framework, naturally fall 
into irreducible representations according to the translational symmetry of 
the system. Thus this method is well adapted to the problem at hand. The 
state-of-the-art of Hartree-Fock (HF) calculations following this approach is 
represented by the CRYSTAL program0. On the other hand, to go beyond 
the HF approach means to allow for electron correlation which is a more or 
less local phenomenon. Here, the use of delocalized Bloch orbitals might not 
be the optimum choice. 

Of course, the latter remark only applies to wavefunction-based meth- 
ods aimed at explicitly including electron correlation at a microscopic level. 
In contrast to that, density-functional theory (DFT) implicitly incorpo- 
rates correlation into a local potential/energy-density, formally still at the 
independent-particle level 0. Very often DFT band-structures are, without 
rigorous theoretical justification, directly identified with the eigenvalue spec- 
trum of the Kohn-Sham (KS) equations. A better theoretical foundation 
has been achieved using the GW approximation which has been developed 
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for calculating quasi-particle excitations in semiconductors and insulators |4]], 
but still the approximations are not fully controlled. 

Wavefunction-based methods for solids suggested in literature so far are 
either of the Quantum- Monte-Carlo or configuration-interaction (CI) type. 
In the latter case, one tries to explicitly exploit the local nature of the corre- 
lations by drawing advantage of the highly developed and accurate ab initio 
quantum-chemical methods|J. Various kinds of such a 'local ansatz' in the 
framework of nonorthogonal orbitals have been proposed by Stollhoff and 
Fulde @ almost two decades ago and yielded a considerable amount of useful 
results since then. More recently the scheme of local increments^] was for- 
mulated using the computationally more attractive orthogonal orbital sets as 
well as standard quantum chemical methods like coupled-cluster (CC). The 
scheme has been applied quite successfully for the determination of ground- 
state properties of a large variety of systems, including elementary^] and 
III-V semiconductors fllPf as well as ionic compounds [|TT| ]. First successful ap- 
plications to correlated band-structure calculations for valence bands have 
also been reported [|12|, |TjJ. However, the main drawback of the incremental 
scheme, both for ground-state properties as well as band-structures, is the 
derivation of the correlation increments from finite-cluster instead of infinite- 
solid calculations. 

By consistently and rigorously taking into account the periodicity of the 
lattice, at all stages of the calculations, we now strive to provide an improved 



procedure. In two recent papers|14|, [15fl , henceforth denoted as I and II, we 



demonstrated the possibility of calculating ground-state properties of ionic 
solids by a procedure which has the novel feature that it directly yields 
localized HF orbitals, as the result of a self-consistent-field (SCF) calculation 
in the linear combination of atomic orbitals (LCAO) framework. In this 
paper, we extend our considerations in a two-fold way. In the first place, we 
focus on the Hartree-Fock ab initio band-structure. Secondly, we enhance the 
numerical stability of our approach, so that calculations of localized orbitals 
become also possible for extended state-of-the-art basis sets approaching the 
HF limit. The localized orbitals at hand, we plan to include correlations in 
future investigations by means of the incremental scheme. 

This paper is organized as follows. In Section 2A we briefly repeat the 
theory within a restricted Hartree-Fock (RHF) framework. A more detailed 
description is given in paper II. Section |2.2| gives the formulae needed to 
obtain the band-structure by our approach. Numerical results for the ground- 
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state energy, the lattice constant, the bulk modulus and the band-structure of 
NaCl are presented in Section |3|. Finally, Section || contains the conclusions. 

2 Theory 

2.1 Hartree-Fock Equations 

Our aim is to directly determine localized (Wannier-type) orbitals for the 
solid. Denoting by |a(Rj) > the Wannier orbitals of a unit cell located at 
lattice vector Rj, the set {|a(Rj) >;a = l,n c ;j = 1,N} spans the occupied 
HF space of the solid. (Here, n c is the number of orbitals per unit cell, 
and N(—> oo) is the total number of unit cells in the solid.) Translational 
symmetry requires 

\a(Ri) >=T(Ri)|a(0) > , (1) 

where T(Rj) is an operator which represents a translation by the lattice vec- 
tor Rj, and |a(0) > denotes the orbitals of an (arbitrarily chosen) reference 
cell with lattice vector R = 0. Assuming now the |a(Rj) > to be orthogonal: 

<«(0)|T(R,)|/3(0) >=5 aP 5 0k , (2) 

the closed-shell HF energy (per unit cell) for the solid can be written as 

n c n c 

E = 2J2 < «(0)|T|a(0) > +2^ < a(0)\U\a{0) > 

a=l a=l 

n c N 

+ £ E( 2 < «(0)/3(R j )|a(0)/3(R i ) > - < a(0) ( 9(R i )|^(R i )a(0) >) 

a,/3=lj=l 

~l" E nuc . (3) 

Here, T and U denote the kinetic-energy and electron-nucleus potential- 
energy operators, respectively, < a(3\a(3 > and < a(3\(3a > represent the 
usual two-electron Coulomb and exchange integrals, and E nuc is the interac- 
tion energy (per unit cell) of the nuclei. The divergent averages of U, E nuc , 
and the Coulomb term are already assumed to have been eliminated from 
Eq. |. 
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Minimization of the |a(Rj) >, under the orthogonality condition, Eq. 0, 
can be shown to lead to the following (modified) Hartree-Fock equations 

(T + U + J E@Jfi-Kf>)+ E £^It(R*) ><7(R*)I)M0) >= e a \a{0) > , 

(3 k£M 7 

(4) 

to be solved in the limit A^ — ► oo. The results are not sensitive to the actual 
values of the A^, as long as they are large enough; in our calculations, we 
chose A^ = 10 4 if. J@ and Kp are the conventional Coulomb and exchange 
operators 

Jp\a> = E i </9(R j )l^l/3(Ri)>l«> \ r ,x 
Kp\a> = Z j </3(R j )\±\a>\/3{R j )> J ' K) 

as discussed in the earlier papers I, II, U, J, K involve infinite lattice sums 
which are rigorously performed, in our treatment, apart from suitable trun- 
cation criteria. For solving Eq. |] self-consistently, we employ a finite-basis-set 
LCAO expansion, with the basis set describing the |a(0) > extending over 
all atoms of a neighborhood M of the reference unit cell. M should consist 
of all unit cells with basis functions having non-negligible overlap to those of 
the reference cell; clearly, the choice of M will also be dictated by the system 
under consideration — the more delocalized the electrons of the system are, 
the larger M will need to be. In our calculations we have typically chosen M 
to include up to third-nearest neighbor unit cells. 



2.2 Band Structure 

After solving the HF equations of the solid for localized Wannier-type orbitals 
|a(Rj) >, we have at our disposal, among other results, the converged Fock 
matrix in the AO basis: 

F p0 , q3 =<p(0)\T+U + ^(2Jp-Kp)\q(^) > . (6) 

(The matrix elements are easily obtained from Eq. |], by subtracting off the 
last (one-electron) term.) 

The band-structure involves matrix elements of the Fock operator in It- 
space. The k-space transform of the AOs \p(Rj) > is 

\ P (k)>=^=Y,^\p(^)>- (7) 

viv R 
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This immediately renders the transformation rule for operators: 

Q w (k) = E^ ikRj (M0)|Q|g(R,)). (8) 

Inserting Q = F (from Eq. |6]) and Q — 1, leads to the k-space Fock and 
overlap matrices, F pq (k) and S pq (k), respectively. A final diagonalization in 
k-space, i.e. solving the generalized eigenvalue equation 

J2F pq (k)C qa (k) = e a (k)Y,S pq (k)C qa (k) (9) 
yields the desired quasiparticle energies for a given k point. 



3 Calculations and results 

We have performed an all-electron calculation for the NaCl crystal adopting 
the experimental face-centered cubic (fee) structure. The unit cell was de- 
fined with sodium at the (0,0,0) position and chlorine at (0,0, a/2), where 
a is the lattice constant. In the following, we will repeatedly compare our 
results to those obtained using the Bloch-orbital-oriented approach of the 
Torino group reported in |T7| for the ground-state properties. For the band- 



structure, we performed calculations using the most recent version CRYS- 
TAL95 of the Torino group 0. In order to establish a strict compari- 



son, we used their most recent basis sets[17], i.e. (15s7p)/[4s3p] for Na and 
(19sllp)/[5s4p] for CI. The p functions were constructed from Gaussian 
lobes [16j which yield an approximate form of Cartesian- type basis functions. 
The number of basis functions per unit cell was 30. To satisfy the orthog- 
onality requirement between wavefunctions located in neighboring cells, a 
neighborhood of 42 unit cells was taken into account according to the pre- 
vious section, which corresponds to up to third-nearest neighbors in the fee 
crystal. Thus the number of basis functions associated with the Wannier 
orbitals is as large as 1290. Note, however, that exploiting equivalence of 
integrals connected by lattice translations effectively reduces the computa- 
tional effort for setting up the Fock matrix so that it is not larger than in a 
corresponding Bloch-orbital-based approach. 

In order to establish a firm basis for a comparison of total energies ob- 
tained within the Bloch-orbital-based approach used in CRYSTAL95 and 
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our present method using Wannier functions we had to guarantee that the 
Coulomb and exchange series in both programs are evaluated with the same 
accuracy. The program CRYSTAL95 offers several input parameters for this 
purpose. In the current work we chose for ITOL1, ITOL3, ITOL4 and ITOL5 
the values of 7, 7, 7 and 15, respectively, and corresponding thresholds in our 
program. This combination of parameters is believed to give well-converged 
results and ensures an accuracy of ~ 1.0 x 10~ 7 atomic units (a.u.) in the 
Fock matrix elements. The resulting accuracy of the total energy is ~ 1 
milliHartree per atom 0. In fact, at the HF equilibrium lattice constant the 
total energies obtained with CRYSTAL95 and our approach agree within one 
mH. A similar good agreement was also observed earlier for other systems 
(cf. I and II). 

For the derivation of the ground-state properties of NaCl, we calculated 
the energy per unit cell for various lattice parameters near the equilibrium 
value (table [I]). We fitted these data points to a cubic polynomial to ob- 
tain the bulk modulus, the equilibrium lattice constant and the equilibrium 
energy. Our results are listed in table |2| together with the corresponding 
numbers obtained by the Torino group with their CRYSTAL program pack- 
age and the experimental values [17|. The lattice energy was calculated with 
respect to the separated Na + and Cl~ ions. To make the comparison with 
the results of the Torino group meaningful, we used the same HF energies 
of -161.67001 a.u. (Na+) and -459.54320 a.u. (CI - ), which can be easily 
reproduced using the basis set reported by them |17|. As can be seen in 
table |2|, the agreement of our results for the lattice constant and the cohe- 
sive energy with the values obtained by the Torino group is excellent. For 
the bulk modulus the disagreement is less than 1 GPa and might at least 
partially be attributed to the use of lobe vs. cartesian basis functions in the 
two approaches, as will be discussed below. 

The band-structure, evaluated at the experimental lattice constant of 
a = 10.526cm, is shown in Fig. [I] for two chosen high symmetry directions in 
the first Brillouin zone. Displayed are the highest three valence bands (the 
upper one being degenerate) as well as the lowest three conduction bands. 
The zero of energy was chosen to coincide with the highest valence point. 
For the sake of comparison, results obtained by using the program package 
CRYSTAL95 of the Torino group with the same basis set are plotted as dotted 
curve together with the results of the present work (solid line). We observe 
an excellent overall agreement when comparing the curves obtained by the 
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two approaches. For the sake of transparency we also give numerical values 
of the quasiparticle energies at the specified symmetry points in table || 

As for the deviations between our results and the ones obtained with 
CRYSTAL95, a remark concerning a technical detail is in order here. As 
can be seen from table |3|, deviations between the two approaches are <2 
mH for the valence bands, whereas for the conduction bands differences of 
up to some ten mH arise for the highest band at the L point, which is well 
outside the overall agreement in the total energy per unit cell given above. 
This can be traced back to the usage of Gaussian lobe basis functions in our 
calculations. To confirm this point we made test calculations on the NaCl 
molecule comparing results with cartesian Gaussian and Gaussian lobe basis 
functions. The comparison of the single-particle energies reveals the same 
features as in the case of the solid. A deviation of some tenth of a mH for 
the highest occupied states is contrasted by a deviation of up to 13 mH in 
the excited states. This helps to clarify the slight deviations in table [| We 
want to point out, however, that the use of Gaussian lobe functions does not 
imply any approximations in our HF scheme, but merely corresponds to the 
use of alternative basis functions. 

4 Conclusions 

An ab initio Hartree-Fock approach for insulating crystals yielding directly 
orbitals in a chemically intuitive localized representation has been applied 
to band-structure calculations as well as the determination of the lattice 
constant, cohesive energy and bulk modulus of sodium chloride NaCl. The 
close agreement between the results obtained using the present approach, 
and corresponding ones obtained using the conventional Bloch-orbital-based 
HF approach, demonstrates that the two schemes are entirely equivalent. 
The advantage of our method is that by taking into account local excitations 
from the Hartree-Fock reference state using conventional quantum-chemical 
methods, one can go beyond the mean-field level and study the influence of 
electron correlations on, e.g., on the binding energy and the band-structure 
of a solid, in an entirely ab initio manner. Work along these lines is currently 
underway in our laboratory. 
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Table 1: Total energies (per unit cell) obtained using our approach for NaCl, 
for different values of lattice constants. Lattice constants are in units of A, 
and energies are in atomic units. 



Lattice Constant 


Total Energy 


5.65 


-621.4930916 


5.693 


-621.4934776 


5.736 


-621.4936939 


5.779 


-621.4937934 


5.822 


-621.4935761 


5.865 


-621.4933206 


5.95 


-621.4929614 
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Table 2: Equilibrium lattice constants (in A), bulk moduli (in GPa), equilib- 
rium energies (in atomic units), and lattice energies (in kcal/mol) for NaCl, 



obtained using our approach and the one reported by the Torino group [[17 
Relevant experimental data are also given for comparison. 



a Ref. 
b Ref. 
c Ref. 



0- 

19]. 

sg. 



Quantity 


Method 


Results 




This Work 


5.785 


Lattice Constant 


Torino 


5.80 




Exp 


5.57 a 




This Work 


23.15 


Bulk Modulus 


Torino 


22.3 




Exp 


28.6 fc 




This work 


-621.494 


Equilibrium Energy 


Torino 


-621.495 




This Work 


176.1 


Lattice Energy 


Torino 


176.1 




Exp 


189.2 C 
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Table 3: Comparison of the quasiparticle energies at high symmetry points for 
the first three occupied and unoccupied bands. Row "a" contains the values 
given by CRYSTAL95, row "b" the corresponding values of this work. The 
energies at the gamma-point are chosen to be zero. All results are in atomic 
units 







1 


2 


3 


4 


5 


6 


L 


a 


-0.094 


-0.014 


-0.014 


0.746 


0.773 


1.096 




b 


-0.095 


-0.014 


-0.014 


0.743 


0.774 


1.047 


r 


a 


0.000 


0.000 


0.000 


0.536 


1.126 


1.126 




b 


0.000 


0.000 


0.000 


0.531 


1.133 


1.133 


X 


a 


-0.084 


-0.031 


-0.031 


0.787 


0.799 


0.974 




b 


-0.086 


-0.032 


-0.032 


0.779 


0.800 


0.978 
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34.0 





Figure 1: Hartree-Fock band structure for NaCl. The highest three valence 
bands (the upper one being degenerate) as well as the lowest three conduction 
bands are displayed. The solid line shows results of the present work. For 
comparison, the results obtained by CRYSTAL95 are shown as dotted lines. 
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